
## if any package is uninstall, use -install.packages("[packagename]") to install first
library(muRL)
library(haven)
library(readr)
library(readxl)
library(rgdal)
library(rgeos)
library(maptools)
library(mapproj)
library(RColorBrewer)
library(ggplot2)
library(usmap)
library(zipcode); data(zipcode)
library(sp); library(sf)
library(albersusa)
shape <- counties_sf("aeqd")

## set local pathname to replication folder
# setwd("[pathname]")

## zips.dta: zipcodes of wealthy respondents
zips <- read_dta("zips.dta")
zips$zip <- ifelse(nchar(as.character(zips$zip)) > 4,
                   as.character(zips$zip),
                   paste0("0", as.character(zips$zip)))

## county-income.csv: average county income from American Community Survey data
acs <- read_csv("county-income.csv")
acs$county <- ifelse(nchar(as.character(acs$region)) > 4,
                     as.character(acs$region),
                     paste0("0", as.character(acs$region)))
names(acs)[names(acs) == "value"] <- "countyincome"

## zip_tract.xlsx: zipcode-tract cross-walk
xwalk <- read_xlsx("zip_tract.xlsx")
xwalk$county <- substring(xwalk$tract, 1, 5)
xwalk <- xwalk[, c("zip", "county")]

shape <- merge(shape, acs, by.x="fips", by.y="county", all.x=TRUE)

## Income categories
cuts <- c(15000, 30000, 40000, 50000, 60000, 750000, 100000, 130000)/1000
shape$qincome <- cut(shape$countyincome/1000,
                     cuts,
                     include.lowest=TRUE)
shape$qincome <- droplevels(shape$qincome)

## geocode_zips.txt: geocoded zipcode data
zipcode <- read.table("geocode_zips.txt",
                      sep="\t", header=TRUE)
zipcode$GEOID <- stringr:::str_pad(zipcode$GEOID, 5, pad="0")
zipcode <- merge(zipcode, zips, by.x="GEOID", by.y="zip")
zipcode <- rbind(zipcode, zipcode[zipcode$GEOID >= 995, ])

# Transform points
spdf <- sf::st_as_sf(zipcode, coords=c("INTPTLONG", "INTPTLAT"))
spdf <- as(spdf, "Spatial")
proj4string(spdf) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
spdf <- points_elided(spdf)

spp <- "+proj=aeqd +lat_0=40.08355112574181 +lon_0=-95.44921875 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"
spdf <- spTransform(spdf, CRS(spp))
shape <- as(shape, "Spatial")

# plot
pdf("FigureA2.pdf")
spplot(shape,
       "qincome",
       col.regions=brewer.pal(9, "OrRd"),
       col="transparent",
       sp.layout=list("sp.points", spdf, pch=1, cex=.8, col="black"),
       ## main="Income distribution by county and location of respondents",
       main="",
       par.settings=list(axis.line=list(col='transparent')),
       colorkey=list(title="Income (in $1,000)"))
dev.off()
